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Preface 

As proven in this work, all orthogonal matrices solve a first order differential equation. The straightfor- 
ward solution of this equation requires n 2 integrations to obtain the elements of the n-th order matrix. There 
are, however, only n(n— 1)/2 independent parameters which determine an orthogonal matrix. The questions 
of choosing them, finding their differential equation and expressing the orthogonal matrix in terms of these 
parameters are considered in this work. Several possibilities which are based on attitude determination in 
three dimensions (3-D) are examined. It is shown that not all 3-D methods have useful extensions to higher 
dimensions. It is also shown why the rate of change of the matrix elements, which are the elements of the 
angular rate vector in 3-D, are the elements of a tensor of the second rank (dyadic) in spaces other than 
three dimensional. It is proven that the 3-D Gibbs vector (or Cayley Parameters) are extendible to other 
dimensions. An algorithm is developed employing the resulting parameters, which are termed Extended Rod- 
rigues Parameters, and numerical results are presented of the application of the algorithm to a fourth order 
matrix. 
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I. INTRODUCTION 


In a recent paper [1] a new algorithm for solving the matrix Riccati equation was introduced. The algo- 
rithm requires the solution of two matrix differential equations. The solution of one of the equations yields a 
diagonal matrix of the eigenvalues of P, the solution matrix of the Riccati equation. The other equation is 

V(t) = W(t)V(t) (1) 

where V is a matrix of the eigenvectors of P. Since P is a real symmetric matrix its eigenvectors are or- 
thonormal, consequently V is an orthonormal matrix. (In the ensuing we will refer to an orthonormal matrix 
as an orthogonal one). The matrix W is a skew-symmetric matrix. (Note that in [1] the order of V and W 
on the right-hand side of (1) is reversed. This difference should cause no difficulty since V is the transpose 
of the corresponding matrix in [1] and W is the negative of its corresponding matrix). 

Let n be the order of the square matrix V. The number of scalar integrations implied by (1) is n 2 ; 
however, the orthogonality of V invokes n(n+l)/2 relations among its elements. Therefore there are really 
only m=n(n-l)/2 independent elements in V. The superfluous computational burden involved in the solution 
of (1) can, then, be reduced by properly defining the m independent parameters of V, solving a differential 
equation only for them and then performing an algebraic computation in order to transform these m elements 
into V. 

We observe that (1) is identical to the famous differential equation of the transformation matrix in the 
three dimensional Euclidean space which is solved on-line for attitude determination of navigation and satel- 
lite systems. That matrix, of course, is also orthogonal, and W is a skew-symmetric matrix whose entries 
are the three components of the angular velocity vector at which the body rotates with respect to some refer- 
ence coordinates. One question that comes immediately to mind is: does (1) always yield a solution which is 
orthogonal? and conversely, do all orthogonal matrices solve such a differential equation? 

The answer to these two questions is formulated in the following two theorems. 

Theorem 1.1: Given equation (1) for t 0 < t < t t where 

WT(t) = -W(t) (2) 


then: 

(I) the matrix V T (t)V(t) is a constant matrix. 

(II) if the initial matrix V(to) is orthogonal, then V(t) is orthogonal too. 
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Proof: 


I [V T (t)V(t)] = V T (t)V(t) + V T (t)V(t) 


( 3 ) 


substituting (1) into (3) yields 


| [V T (t)V(t)] = V T (t)W T (t)V(t) + V T (t)W(t) V(t) 
and when (2) is substituted into (4), it is seen that 

| [V T (t)V(t)] = 0 


( 4 ) 


( 5 ) 


Consequently 


V T (t)V(t) = Const. 


(6) 


and thus (I) has been proven. 

Now when V(t 0 ) is orthogonal, then 


V T (t 0 )V(t 0 ) = I 

(where I denotes the identity matrix) and due to (6) also 

V T (t)V(t) = I 


which proves assertion (II). ■ 

Theorem 1.2: Any time varying orthogonal matrix, V(t), satisfies the matrix differential equation 

V(t) = W(t)V(t) (7) 

where 

W T (t) = -W(t) (8) 


Proof: Since V(t) is orthogonal 


V(t) = V(t)V T (t)V(t) 


(9) 


Denote 


W(t) = V (t)V T (t) 
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then (9) can be written as 


V(t) = W(t)V(t) 


( 10 ) 


which is (7). 

Using (10) we write 

V T (t)V(t) + V T (t)V(t) = V T (t)W T (t)V(t) + V T (t)W(t)V(t) = V T (t) [WT(t) + W(t)] V(t) (11) 
The left-hand side of (10) is the time derivative of V T (t)V(t) hence (10) can be written as 

| [V T (t)V(t)] = V T (t) [WT(t) + W(t)] V(t) (12) 


But 


V T (t)V(t) = I 

hence the left-hand side of (11) is zero which implies that 

WT(t) = -W(t) 


as stated in (8). This completes the proof. H 

In view of the preceding, it is realized that the problem we are concerned with is an extension of the 
three dimensional attitude determination problem and conversely, the latter is a special case of the problem 
at hand. It is interesting to investigate the correspondence of the various elements involved in three dimen- 
sional attitude determination with the eventual solution and features of our present problem. For this reason 
the pertinent background material of attitude determination will be reviewed in Section HI following a for- 
mal definition of the problem in the next section. In Section IV we discuss a possible solution using Extend- 
ed Euler Angles followed, in Section V, by an introduction of the chosen Extended Rodrigues Parameters 
solution. In Section VI we probe the issue of presenting angular rate in n-D and in Section VII we discuss 
numerical issues involved in the implementation of the solution. Numerical results are then presented and 
conclusions are drawn in Section VIII. 


II. PROBLEM STATEMENT 

We state our problem as follows. Given the matrix differential equation 

V(t) = W(t)V(t) 

in which W is a skew symmetric matrix and for which the initial matrix V(t 0 ) is known to be orthogonal, 
find the following: 

a) m=n(n— 1)/2 parameters which unambiguously define V, 
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b) the differential equation needed to be solved in order to compute these parameters, 

c) the functional relations between the parameters and V which will enable the computation of V based 
on the parameters, and 

d) a simple algorithm to implement the solution of the differential equation as well as the computation 
of V. 


in. BACKGROUND IN THREE DIMENSIONAL SPACE 


Euler Angles 

The best known parameters describing a 3-D rotation and the resulting transformation matrix are Euler 
Angles [2-6]. Three such angles are necessary and sufficient to describe any transformation from one Car- 
tesian coordinate system to any other one. There are 12 sequences of 3 right-hand Euler Angle rotation 
sequences. If for example one chooses the sequence z-y-x rotations by the respective angles p, t and f, then 
the corresponding differential equations of the Euler Angles are: 


p = (w y sin f + w z cos f)/cos t 

(13. a) 

t = w y cos f — w z sin f 

(13.b) 

f = w x + tan t (w y sin f + w z cos f) 

(13.c) 


where w x , w y and w z are the three components of the angular rate vector at which the final coordinate sys- 
tem turns with respect to the initial one when this vector is resolved in the final system. The transformation 
matrix, D, which transforms vectors from the initial coordinate system into the rotated one is computable 
using the solution of (13) in the following expression 



cp ct 

sp ct 

—st 

D =• 

— sp cf -fcp St sf 

cp cf +sp st sf 

ct sf 


sp sf +sp St cf 

— cp sf +sp St cf 

ct cf 


where s denotes the sine and c denotes the cosine function. 

We note two shortcomings of this method. First, we run into a singularity problem as t approaches 90° 
or —90° and, secondly we need to compute trigonometric functions. For this reason the use of Quaternions 
is usually preferred. 

Quaternion 

Quaternions consist of 4 elements; that is, the Quaternion is a 4 parameter rotation specifier [2,5,6]. 
One parameter is, of course, superfluous but this is acceptable since, using Quaternions, the two aforemen- 
tioned shortcomings, involved in the usage of Euler Angles, are eliminated. Denote the 4 elements of the 
Quaternion of rotation by q 0 , qi, q 2 and q 3 then the differential equation of the Quaternion elements is 
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“qo“ 


“0 

-W x 

-w y 

-w z 


"qo~ 

d 

qi 

_ i 

W x 

0 

W 2 

-Wy 


qi 

dt 

q2 

2 

w y 

-w z 

0 

Wx 


q2 


_q 3 _ 


_W Z 

Wy 

-W x 

0 _ 


-q 3 - 


( 14 ) 


The solution of (14) yields the components of the Quaternion which can be used to compute D as follows 


qo 2 + qi 2 - q 2 2 - q3 2 2 (qiq 2 - q 0 q 3 ) 


D = 


2(qiq2 + qoq 3 ) 
2(q 3 qi ~ qoq 2 ) 


qo 2 - qi 2 + q 2 2 - q 3 2 
2(q 2 q 3 + q 0 qi) 


2(q 3 qi + qoq 2 ) 

2(q 2 q 3 - qoqi) 
qo 2 - qi 2 - q 2 2 + q 3 2 _ 


The Quaternion of rotation is based on Euler’s theorem which states that any orientation of a 3-D Cartesian 
coordinate system with respect to any reference system can be obtained by a single rotation of the initial 
coordinate system about an axis fixed in both systems. Let the positive direction (according to the right-hand 

A 

rule) of this axis be denoted by a unit vector f and the rotation angle by f, then the components and the 
magnitude of the rotation vector ff (also known as Euler Vector) are used to define the Quaternion as 
follows 


qo = cos(f/2) ; qi = sin(f/2)f x /f 


q 2 = sin(f/2)f y /f ; q 3 = sin(f/2)f z /f 

where f,, i=x,y,z, are the 3 components of the rotation vector. The Quaternion is, then, a 4 component ele- 
ment constructed on a 3 component vector. 

Rodrigues Parameters 


Another 3 parameter representation of 3-D rotations is due to Rodrigues [7]. Denote the parameters by 
gi, g 2 and g 3 , then the differential equation which these parameters satisfy is [7,8,6] 



— “ 


d 

gi 

1 

dt 

g2 

2 


g 3 



1 + gi 2 
“g 3 + glg2 
g2 + glg 3 


g 3 + glg2 
1 + g2 2 
“gl + g2g 3 


The solution of (15) can, then, be used to compute D as follows 


“ 


— 

g2 + glg 3 


W x 

gl + g2g 3 


Wy 

1 + g 3 2 


W z 


(15) 


d = 1 + gi 2 + g 2 2 + g 3 2 



1 + 


gl 2 ~ g2 2 - g 3 3 


2(glg2 - g 3 ) 

- gi 2 + g2 2 - g 3 2 
2(g2g 3 + gl) 


2(glg 3 + g2) 
2(g 2 g 3 “ gi) 

- gl 2 - g2 2 + g 3 2 


2(gig 2 + g 3 ) 
2(gig 3 - g 2 ) 
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The relationship between the Rodrigues Parameters and the rotation vector are 

gt = tan(f/2)f x /f ; g 2 = tan(f/2)f y /f ; g 3 = tan(f/2)f z /f 

Since both the Quaternion of rotation and Rodrigues Parameters are based in a similar manner on the rota- 
tion vector, there is a rather simple relationship between them; namely, gi = qj/q 0 for i= 1,2,3. 

The preceding equations for the time change of the Rodrigues Parameters and for converting the 
parameters into D can be cast in matrix form as follows [9,10]. Define a G matrix such that 


and, similarly a W matrix 


then 



(16) 


(17) 


G = - |(I + G)W(I - G) (18) 

D = (I - G) (I + G)- 1 (19) 

where I is the identity matrix. Like with the 3 parameter Euler Angle representation, here too, singularity 
may occur whenever the size of the rotation vector reaches a magnitude of 180°. 

After having discussed the possible solutions to the problem in 3-D we will consider, next, the possibili- 
ty of extending these solutions to n-D (whenever mentioning n dimensional spaces we mean Euclidean 
spaces whose dimension n:£3). 

IV. POSSIBLE SOLUTIONS 

Extended Euler Angles 

When trying to solve our problem (as defined in Section II) the first question that comes to mind is: can 
the Euler Angle parametrization presented in the preceding section be extended to higher dimensional Eu- 
clidean spaces? As it turns out [11], Euler himself showed that this was possible. This was also shown later 
by Lagrange [12]. (See also Jacobi’s observation on their and other’s work [13]). However, the use of the 
Extended Euler Angles for n > 3 is cumbersome since, for calculating V, the sine and cosine functions of 
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m=n(n— 1)/2 angles must be computed, these functions have to be multiplied through in a long string of 
multiplications, and the resultant products have to be added and subtracted. For n=4, for example, the 1,1 
element of Y is 


vi i = cos ai cos a 3 cos as + sin sin a 4 sin as 

and there are 16 elements, all equally long, in V. When compared with the simplicity of the solution which 
we will eventually choose the complexity of the present one will be striking. Moreover, to complete the al- 
gorithm it is necessary to find the differential equations governing the Extended Euler Angles and solve 
them. Merely finding the equations, let alone solving them, is a formidable task. As an example for the 
work involved in deriving those equations, consider the following approach. Let the Extended Euler Angles 

be denoted by ai, a 2 , , a m . Denote the column vector whose elements are these angles by a. We 

may express V as a product of the individual matrices V(a ; ) of the transformation matrix related to a single 
angle a,, thus 


m 

v (a) = V(aj) (20) 

i = l 

Differentiation of (20) yields 


V(a) = £ a j 
j=l 

On the other hand (20) and (1) yield 


V(a) = W n V(a,) (22) 

i = l 

equating the right-hand sides of (21) and (22) yields m equations in aj. After cumbersome manipulations we 
obtain the required m differential equations for aj, j = l,2,...,m whose solution yields a, the elements 
of which are needed in order to compute V. We conclude that finding the differential"equations for the 
Extended Euler Angles, solving them, and then using the solutions to compute the corresponding V matrix, 
while possible, is indeed a formidable task which we reject in favor of the method which we will eventually 
select. 

Extended Quaternion 

The use of the quaternion of rotation in 3-D is motivated by the following considerations. It does not 
suffer from singularities, it does not require the computation of trigonometric functions, it has a simple 
linear differential equation and a simple geometric interpretation related to the rotation vector. Finally, the 
only price paid for using it, is the need to deal with 4 (rather than 3) parameters. Because of these merits, 
one is motivated to try to extend the notion of quaternions to n-D. This approach though does not seem to 
yield a non-singular parametrization even if one is willing to use m+1 parameters to define an extended 
quaternion. 
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Of the three 3-D parametrization methods reviewed in Section HI only the Rodrigues Parameters are ex- 
tendible to a compact easily implementable algorithm. This will be shown in the next section. 

V. EXTENDED RODRIGUES PARAMETERS 

We start the presentation of this parametrization method in n-D with two lemmas which will be helpful 
in the ensuing. 

Lemma V.l: Let A be an nxn matrix, then the matrix (I + A) is invertible if none of the eigenvalues of A 
is equal to —1. 


Proof: 

The eigenvalues , i = l,2,...,n of (I + A) are the roots of the polynomial 

| (I + A) - bl | =0 


which can be written as 


A - (b - 1)1 1 = 0 


(23) 


(24) 


or 


where 


| A - al | = 0 


a = b - 1 


(25) 


(26) 


The condition for (I + A) to be invertible is b, 4 0, i = l,2,...,n or, in view of (26), a; 4 — 1, 
i = 1,2, ...,n., But in view of (25), aj are the eigenvalues of A. This ends the proof. 

Lemma V.2: Let (1 + A) -1 exist and let 


B = (I - A) (I + A)- 1 


then (1 + B) is invertible. 


Proof: 


(I + B) = I + (I -A) (I + A)" 1 

= (I + A) (I + A)" 1 + (I - A) (I + A)” 1 
= 2(1 + A)- 1 
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Obviously, (I + A) -1 has an inverse which is (I + A), thus 

(I + B)- 1 = (I + A)/2. ■ 

With these lemmas on hand we can proceed and prove the following theorem. 

Theorem V.l: Let V be an n-th order orthogonal matrix with none of its eigenvalues equal to —1; then 

I) there exists a matrix G defined as follows 

G = (I -V) (I + V)" 1 (27) 


II) G is skew-symmetric 

III) V is the following function of G 


V = (I - G) (I + G)" 1 (28) 

IV) The rate of change of G is given by 

G = - |(I + G) W (I + G) T (29) 

where W is the matrix defined by (1). 

Proof: From lemma V.l the matrix (I + V) has an inverse; thus G as defined in (27) exists. To show that 
G is skew-symmetric use (27) to write 


G T = (I + V) -T (I - V) T 

where — T is the inverse of the transpose (or vice-versa). Using the last equation and the orthogonality of V 
we observe that 

G T = (I + V T )- 1 (I - V T ) = (V T V + V T )-‘(V T V - V T ) 

= [V T (V + I)]-tV T (V - I) 

= (I + V)-1VV T (V - I) = (I + V)-!(V - I) 

= -(I + V)-![2I - a + V)] = - 2(1 + V) -1 + I (30) 

Now 

-2(1 + V)" 1 + I = -2a + V)" 1 + (1+ V) (I + V)" 1 
= [-21 + (I + V)] (I + V)" 1 

= - (I - V) (I + V)- 1 = -G (31) 


9 



Substitution of (31) into (30) yields the result G T = -G, i.e. G is skew-symmetric. 

From lemma V.2, (I + G) is invertible which gives legitimacy to the right-hand side of (28). To prove 
the truth of (28) re-write (30) as 

G T = I - 2(1 + V)' 1 

hence 

G = I - 2(1 + V 1 )" 1 

from which we obtain 

I - G = 2(1 + V 1 )- 1 (32) 

and 

I + G = 21 - 2(1 + V 7 )" 1 (33) 

We can further write 

21 - 2(1 + V 7 )" 1 = 2(1 + V 7 ) (I + V 7 )- 1 - 2(1 + V 7 )" 1 
= 2V T (I + V 7 )- 1 

thus 

I + G = 2V T (I + V 7 )" 1 
and 

(I + G)- 1 = |(I + V 7 )V (34) 

Substitution of (32) and (34) in the right-hand side of (28) yields the proof of III. 

To prove (29) differentiate (27) 

G = — V(I + V)" 1 - (I - V) (I + V) -1 V(I + V)- 1 
= -[I + (I - V) (I + V)- l ]V(I + V)- 1 
Substitute (27) in the last equation to obtain 

G = —(I + G)V(I + V)- 1 
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Using (1) the last equation can be written as 


G = -(I + G)WV(I + V)" 1 
Substitution of (28) into the last equation yields 

G = -(I + G)W(I - G) (I + G) -1 [I + (I - G) (I + G)" 1 ]- 1 (35) 

The expression in the brackets can be written as follows 

I + (I - G) (I + G)- 1 = (I + G) (I + G)" 1 + (I - G) (I + G)- 1 = 2(1 + G) -1 
therefore (35) can be written as 

G = -(I + G)W(I -G) (I + G) _1 [2(I + G)- 1 ]' 1 = - |(I + G)W(I - G) 


and since G is skew-symmetric the last equation can be written also as 

G = - 1(1 + G)W(I + G) T 


which ends the proof. H 

Note, from lemma V.l, that the condition for the invertibility of (I + G) is that it has no eigenvalues at 
-1, which is analogous to the condition for (I + V) to be invertible, i.e. that V has no eigenvalues at -1. 
However while V always exists, G does not exist when V has an eigenvalue at —1. The parametrization of 
V by G fails when the latter is the case. However this can be overcome as will be shown in Section VII. 

The parametrization of V by the Extended Rodrigues Parameters is n-dimensional since the foregoing 
proofs were not restricted to any value of n, nor did they hinge on a rotation vector or any other geometric 
quality in n-D. In fact, the Extended Rodrigues Parameters, which are the elements of G, are the answer to 
the first three parts of our problem as posed in Section II. That is, we found m parameters which define the 
n-dimensional orthogonal matrix, V. We also found a first order differential equation for G, and we showed 
how to calculate V, once G is found. 

What is needed to fully answer our problem is a simple algorithm to implement the solution; this will 
be presented in Section VII. For now, after having obtained a parametrization in n-D, we are prepared to 
discuss the meaning of the skew-symmetric matrix, W, its geometric interpretation, and the difference be- 
tween W in 3 and in n-D. 


11 



VI. ANGULAR RATE IN n-D 


Recall (1) 


V(t) = W(t)V(t) (1) 

The matrix V can be viewed as a transformation matrix which transforms vector components in an n-D Eu- 
clidean space. In particular it transforms a set of unit vectors, which form a Cartesian coordinate system, to 
another such set. Let us denote the former as the initial coordinate system and the latter as the final one. 

The rows of V are components of unit vectors of the initial set resolved in the final Cartesian coordinate 
system such that Vj j is the i-th component in the final system of the j-th unit vector of the initial coordinate 
system. From (1) 


n 



k=l 


hence Wj k is the relative weight that the k-th component in the final system, of a unit vector in the initial 
system, has on the rate of change of the i-th component in the final system of the same unit vector in the 
initial system. Note that this weight is independent of j; i.e. of which unit vector in the initial system we 
consider. To give wy a more descriptive interpretation and to see the role of W more clearly, consider the 
3-D case where, for example 


*3,1 = w 3>1 v u + w 3>2 v 2il (36) 

(note that the term w 3 3 v 31 was dropped since w 3 3 = 0 for skew-symmetric W). In 3-D (36) can be written as 

*3,1 = W 2 v u - W!V 2)1 (37) 

where W[ and w 2 are the respective angular rates at which the final coordinate system instantaneously rotates 
about its 1 and 2 axes. The components w i5 i=l,2,3, are those of the 3-D angular rate vector describing the 
instantaneous rotation of the final system. In 3-D w; is also the angular rate at which the j axis turns 
towards the k axis, and so on in a cyclic manner for wj and w k . Indeed a comparison between (36) and (37) 
reveals that 


w 2 — W 3>1 
— Wj = W 3>2 

We conclude that the following can be said about W in 3-D 

(A) The elements of W are angular rates. 

(B) Each components of W is a rate of turn of one coordinate axis towards another 
such that w p q is the angular rate at which axis p turns towards axis q. Obviously, 


12 



(C) Both the p and the q axes turn at the angular rate w p q about the third axis r. 

(D) The elements of W are components of an angular rate vector. 

When we turn now to n-D, we realize that the preceding observation cannot be fully extended from 3 to 
n-D. In n-D, W has m = n(n - l)/2 independent components such that the elements of W cannot be com- 
ponents of a rate vector whose number is necessarily only n. We cannot, therefore, consider the elements of 
W as angular rates about (coordinate) axes. Consequently, of the four features of the elements of W in 3-D, 
mentioned above, the only ones which also prevail in n-D are (A) and (B). 

Realizing that the angular rates in n-D cannot be described by a vector, one is motivated to examine the 
possibility of expressing the angular rate by a tensor. To accomplish that, choose one, say the i-th, column 

of V(t) and the i-th column of V(t) and denote them correspondingly by v and v such that v - [v t , v 2 , 

v n ] T and v = [vi, v 2 , , v,J T . Using their components express them as vectors in the same arbitrarily 

chosen coordinate system such that 


v = i^i + i 2 v 2 + + i n v n 

where ii, i 2 , ..., i n are unit vectors along the coordinate axes 1, 2, ..., n respectively. Similarly 

v = ijV! + i 2 v 2 + + i n v n 

Define a tensor of the second rank, W, using the elements of W as follows 

W = ijiiO + iii 2 Wi >2 + + Mn w i,n 

+ i 2 iiw 2 j + i 2 i 2 0 + + i2in w 2,n 


+ inilWn.l + i n i2W n ,2 + + WnO 


then obviously 


v = W v 

that is, when the angular rate components are treated as elements of a tensor of the second rank, (1) is fully 
satisfied. A tensor of the second rank is also known as dyadic [14], 

The fact that the angular rate in 3-D is basically a tensor is known [8,15] but is not reflected in the 
applied literature. The reason for it stems, perhaps, from the unique possibility of expressing angular rates 
in 3-D by a vector such that its description as a tensor might have been perceived merely as a philosophical 
formalism. (Even when treated as a tensor, the angular rate is usually that of a 3-D coordinate system). In- 
deed, the creation, in 3-D, of the so called “vector cross-product matrix” based on the angular velocity 
vector is conceived as a useful gimmick rather than a restoration of the true mathematical description of the 
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angular rate. So far, the consideration of angular rates in dimensions higher than 3 probably was not re- 
quired nor known. Thus it was not recognized that in higher dimensions the angular rate cannot be described 
by a vector but must be described as another entity, and that the ability to describe it in 3-D by a vector is 
just a matter of good fortune. (In fact, even in 2-D the angular rate is not truly expressible as a vector. This 
is evident when we note that the expression of rotation in a plane by a vector normal to it is necessarily a 
3-D expression. The correct and only 2-D expression is 

v = ( ift 2 w - i 2 ijw )v 




— * — 


— — 

Vl 


0 w 


Vl 

v 2 


-w 0 


v 2 

— — 


— __ 




where the first expression is in a tensor form and the second is in a matrix form). Another possible cause 
for the disregard of the fact that angular rate is a tensor stems from the fact that the tensor of the second 
rank, that is, the dyadic, is replaceable by a matrix (as demonstrated in the last 2-D representation and in 
equation 1). Therefore all practical work in any dimension can be carried out without resorting to the tensor 
concept. 


After having cleared the issue of angular rate representation we are prepared to consider the implemen- 
tation of the algorithm for solving (1) using the Extended Rodrigues Parameters, thereby solving our 
problem in its entirety. 


VII. NUMERICAL IMPLEMENTATION 

Recall the differential equation (1) 


V = WV ( 

in which W is given. We wish to solve (1) using the extended Rodrigues Parameters. The solution process 
requires first the solution of 


G = - |(I + G)W(I + G) T (29) 

and then the computation of V according to (28) 

V = (I - G) (I + G)-i (28) 

There are two caveats which we have to be alerted to. One of them is the non-existence of G when V has 
an eigenvalue at ~ 1, and the other is the need to invert the matrix (I 4- G), which may be so burdensome 
as to render the whole approach inefficient in comparison with the direct solution of (1). The first problem 
can be easily avoided if we can keep the elements of G small, for then, as can be readily seen from (28), V 
is close to I whose eigenvalues are all equal to +1. That is, if we are free to control its size, we can always 
choose G so small as to make the eigenvalues of V as close to +1 (and thus as far from -I) as we wish. 
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Indeed, we are able to control the magnitude of G. The ability to do it is based on the following propo- 
sition. 


Proposition: Given the differential equation of (1) 

V(t) = W(t)V(t) (1) 

with the initial condition V(to) where V(t 0 ) is orthogonal, then V(t), the solution of 
(1) at time t > t 0 , can be written as a product of two matrices as follows 

V(t) = V(t,t 0 )V(t 0 ) (38) 

where V(t,t 0 ) is the solution of (1) at time t given the initial condition V(t 0 ,t 0 ) = I. 

Proof: Since V(to) is orthogonal it always has an inverse. Therefore one can always compute a matrix 

V(t,to) = V(t)V T (to) (39) 

such that (38) holds. Now if (38) is differentiated with respect to time the following is obtained. 

V(t) = V(t,t 0 )V(t 0 ) 

Equating the right-hand side of the last equation to that of (1) and using (38) results in 

V(t,to)V(t 0 ) = W(t)V(t,t 0 )V(t o ) 

Since V(to) is invertible, the last equation yields ■ 

V(t,to) = W(t)V(t,t 0 ) (40) 

hence V(t,t 0 ) solves (1). Finally setting t in (39) to t 0 results in 

V(to,t 0 ) = I 


which ends the proof 

In computing V(t) we make use of the last proposition and instead of computing V(t) directly we com- 
pute V(t,t 0 ) from time t 0 to t and then use (40) to compute V(t). Actually instead of computing V(t,to) we 
use (29) to compute G, the parametrization of V(t,to), from time t 0 to t with the initial condition G(t 0 ) = 0 
which corresponds to V(t,to) = I. The computation of G is stopped periodically at, say, ti and V(ti) is com- 
puted according to (28) yielding 


Vfa.t,,) = [I - G(t,)][I + GOO]’ 1 
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and then V(ti) is computed using (38) as follows 


veto = vo^yvdo) 

Next the computation of V(t) proceeds into the following time interval using the same algorithm that 
produced V(ti) once V(t 0 ) was given. We start, of course, with the initial condition G(ti) = 0 which cor- 
responds to V(t,ti) = I. Using this algorithm we proceed to compute G and V at times t 2 , t 3 t k . By 

properly choosing the size of the intervals t 2 — ti, t 3 — t 2 , t k — t k _i we can impose an upper bound 
on G which can practically be as small as we wish. We term the operation of resetting the value of V and G 
at the beginning of an interval reset operation. 

The foregoing policy rids us of the singularity problem. In fact, if singularity were the only issue, one 
can choose the time intervals t ; — ti_i quite large and still not encounter singularity. However, we are still 
left with the second problem mentioned before; namely, the inversion of [I + G(tj)]. We overcome this 
problem by approximately the inverse without really performing any matrix inversion. Before discussing the 
options for approximating this inverse we list without proof two well known theorems (e.g. Ref. 16 p. 129) 
needed in the ensuing. 

00 

Theorem VII. 1: Let G be a square matrix then the series (— l^G* converges to 

i=o 

(I + G) -1 if all the eigenvalue of G lie inside the unit circle about the ori- 
gin of the complex plane. 

Theorem VII.2: Denote the elements of the nxn matrix G by gy. If the sums 


I Si J I 
j = l 

are all less than 1, or if 


i = 1, 2, ..., n 


n 



i = l 


are all less than 1, then all the eigenvalues of G lie inside the unit circle 
about the origin of the complex plane. 

The algorithm we use to approximate the inverse of (I — G) is based on the fact that if the matrix Xj is 
a good approximation of the inverse of some matrix A then a better approximation, X; + i, can be obtained 
using the Newton-Raphson-type iteration [16 p. 52, 17] 


Xi + ! = Xi (21 - AXi) (41) 

This algorithm converges if and only if the eigenvalues of I — AXj are all of absolute value less than 1 [16, 
p. 52]. If indeed X, is almost the inverse of A then this condition is met. If now V is computed without 
reset taking place at the end of the previous time increment, then we use as a first approximation of 
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[I + G(t)] ~ 1 , the value used as an inverse at the previous time point. This is based on the presumption that 
the time increments of the integration are small enough such that the change of the inverse is small too, 
hence its previous accurate value can serve now as an approximate value. If, however, reset did take place 
at the previous time point then G was set to zero and the previous inverse of I + G is simply I. For the 
sake of computation reduction it is desired to keep at minimum the number of iterations used to compute an 
accurate inverse. Normally one iteration is sufficient. However, when reset takes place and consequently the 
previous inverse of I + G (i.e. the inverse of I + G 0 for G 0 = 0) is taken as I then a single iteration 
produces 


(I + GO-i ~ I - Gj (42) 

where Gj is G at the present time. If, however, we enter the iteration with the value X 0 = I - Gi then, 
due to the quadratic convergence characteristic of the process, a single iteration produces 

(I + GO -1 ~ I - Gi + Gj 2 - Gj 3 (43) 

obviously the approximate inverse given in (43) is more accurate than that of (42) since it contains more 
terms of the series which expresses the inverse of (I + GO- Note that the series generated by (41) con- 
verges since due to the reset operation, G is kept at a very small value such that the condition of theorem 
VII. 2 is met. Thus the eigenvalues of G are in the unit circle which, in view of theorem VII. 1, assures 
convergence. 

Another point of interest is the ability to use an alternate equation for computing G. From (31) it is ob- 
vious that 


2(1 + V)- 1 = I + G 


which yields 


V = 2(1 + G)-i - I (44) 

The computation of V using (44) is simpler than when (28) is used. However, if the reset operation took 
place at the previous time point then the use of (28) at the present time point yields better results. This is 
evident in particular when the previous inverse, (I + Gj _ 0 -1 , is approximated by I. In this case the use of 
(41) yields the approximation of (42) for which the use of (44) yields 


V(Mi _ j) = I - 2Gj 


whereas the use of (28) yields 


V(ti,tj -0 = 1- 2Gj + Gi 3 

which is more accurate than the preceding result. Even when the approximation of (43) is used, the use of 
(28) yields better results than that obtained using (44). Then, however, the difference is smaller since the 
term of the series which is being added is smaller than the added term in the previous case which was Gi 2 . 
If, of course, an exact inverse is used then the use of (44) rather than (28) is preferable since then the com- 
putation of V(t,, t; _ i) is simplified without the penalty of accuracy degradation. 
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The algorithm which results from the preceding considerations is shown in Table I. Note that (28) rather 
than (44) is implemented for the reasons discussed above. 

Table I 

Given: V(t 0 ) = V 0 and W(t) 

(1) Set the initial conditions G(t 0 ) = 0 and i = l. Consider t 0 as a reset point. 

(2) Solve G(t) = -|[I + G(t)]W(t)[I + G(t)] T from t^ to q. 

(3) Initialize X, as follows: 

(a) If reset took place at the end of the preceding cycle, compute X, = I — G(tj). 

(b) If reset didn’t take place at the preceding cycle compute X; = X*; _ i 

(4) Compute A; = I + G(t ; ) and Xi* = X; (21 - AjXj) 

(5) If reset is requested or if current time is equal to the final time perform a reset as follows (other- 
wise go to 6): 

(a) Vft.ti-!) = [I - G(ti)]Xi* 

(b) V( ti ) = 

(6) If the current time is equal to the final time go to step (8). 

(7) If reset was performed in this cycle go to (1). Otherwise increase i by 1 and go to (2). 

(8) Stop. 


If one chooses to perform reset after each integration step of G then the computation of V(ti,t 0 ) as given 
in Table I produces 


V(ti,t 0 ) = I - 2Gj + 2Gj 2 - 2G ; 3 + Gj 4 (45) 

where Gj = G(tj). This is a truncated series of the expression for V as a function of G given in (28) with 
the special feature that the last term in the series lacks the multiplier 2. A more computationally efficient 
algorithm than that is 
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V(ti,t 0 ) = I - G;{21 - Gj[2I - Gj(2I - Gi)]} 


(46) 


or better yet, if Gj 4 is added to (45) to generate the true truncation of the series expansion of (I — G,) 
(I + G;) _1 then the new expression can be written as 


V(ti,t 0 ) = I - 2Gi{I - Gj[I - Gi(I - G;)]} (47) 

Consequently one can use either the algorithm of Table I as is or compute V(ti,t 0 ) using either (45) or (46) 
or (47). These, however are not the only possible variants of the algorithm. As a result of the discussions 
presented in this section it is clear that one has the following additional choices: 

• Perform or not perform resets. 

• Use more terms of the series 

00 

V(ti,ti- ,) = I + 2^(“l) n Gi n 

n = l 


• Use either (28) or (44) to compute V(ti,t 0 ) from Gj. 
The choices should correspond to the particular problem on hand. 

As an example we ran a 4th dimensional case where 


V(0) = I ; W(t) = 


0 

-0.1 

-1.0 

-7.5 

0.1 

0 

3.0 

0 

1.0 

-3.0 

0 

-0.9 

7.5 

0 

0.9 

0 


sin(6.28t) 


the initial time to = 0. 
the final time t f = lsec 
the integration interval dt = O.OOlsec 

The algorithm used in the solution of V was the one given in Table I where reset was performed after each 
integration step. Equation (1) was solved to yield a reference with which the algorithm output was com- 
pared. The reference matrix was denoted by V r and the one generated by the algorithm was denoted by V. 
The integration routine which was used to solve the differential equation for V r as well as for G was a 4-th 
order Runge-Kutta routine. The difference matrix between the two solutions was computed and denoted by 
E = V — V r . A scalar which constitutes a measure of the size of the error was defined as follows 

e = [Tr{EE T }] 1/2 
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The scalar e is the square root of the sum of the squares of the elements of E. The results at t = 0.5sec 
were: 


V r 


— .72765515E+00 
.10217642E— 01 

— .13935294E+00 
.67156112E+00 


— .72765512E+00 
.10217638E— 01 

— .13935294E+00 
.671561 10E+00 


.15285696E+00 

.58373643E+00 

— .79737729E+00 

— .87171959E— 02 


.15285696E+00 

.58373643E+00 

— .79737729E+00 

— .87171923E— 02 


— .24387237E+00 
.79194147E+00 
.53481405E+00 

— .16531458E+00 


— .24387236E+00 
.79194147E+00 
.5348 HOSE +00 

— .16531458E+00 


— . 62263874E + 00 

— .17881859E+00 

— .24237192E+00 

— .72221933E+00 


— .62263872E+00 

— .17881859E+00 

— .24237191E+00 
-.72221930E+00 


E 


.29723949E— 07 

— .35405291E— 08 
.42739559E— 09 

— .25932268E— 07 


— .18713478E— 09 

— .25268088E— 09 
.12413215E— 08 
.35672249E— 08 


.83326148E— 08 
— .11170641E— 08 
.86561824E— 09 
.95984923E— 09 


.24813968E— 07 
.71958844E-09 
.83593105E— 08 
.29599694E— 07 


e = .56724776E— 07 

As mentioned earlier the algorithm of Table I with a reset at the end of each integration cycle amounts 
to the use of (46) in the computation of V(ti,to). As suggested, (47) can be used instead. In Table II we 
show a comparison between the use of (46) and (47) for different series lengths. The table presents the error 
measure, e, for the two series truncated after different powers, n, of G. The error measure was recorded at 
t = 0.5sec for at that point, which is half the period of the oscillating W, the value of e is the highest dur- 
ing the first period, i.e., in the domain 0. < t < l.sec. As can be seen from the results, algorithm 2 is su- 
perior. It can be also seen that there is a distinct power beyond which the addition of more terms yields 
little return. In view of these conclusions we recommend the use of the algorithm listed in Table III which 
in fact was used in the first example. 
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Table II 



n 

1 

2 

3 

4 

5 

1. V(ti,to) = I - 2Gi + 2Gj 2 - ... Gi" 

.17 E 01 

.52 E— 02 

.17 E— 04 

.57 E— 07 

.13 E-09 

2. V(ti,t 0 ) = I - 2Gj + 2Gj 2 - .. 2Gi" 

.10 E— 01 

.34 E— 04 

.11 E— 06 

.33 E— 09 

.63 E— 10 


Table III 

Given: V(to) = V 0 and W(t) 

(1) Initialize i = 0 

(2) Set the initial condition G(ti) = 0. 

(3) Solve G(t) = -|[I + G(t)]W(t)[I+ G(t)] T from h to t; + j. 

(4) Compute 

V(ti + !,tj) = I — 2G; + j{I - Gj + iP - Gi + i(I - Gi + i)]} 
v(ti + 1) = V(ti + Ltovcti) 

(5) If the current time is smaller than the final time go back to step (2) and increase all indices by 1 , 
otherwise STOP. 
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vm. CONCLUSIONS 


This work addressed the problem of solving the first order differential equation, which every orthogonal 
matrix satisfies, using the minimum number of parameters necessary to uniquely determine the matrix. The 
major question was: which are the parameters that determine this matrix. The other questions were: what is 
the differential equation which one has to solve in order to find the parameters, and: once the parameters 
are found, how to use them in order to find the corresponding matrix. All these questions were answered 
and several algorithms for computing the orthogonal matrix via the parameters were suggested and inves- 
tigated. 

In search for solutions the familiar special 3-D case was examined with the purpose of extending the 
methods used there to the general n-D case. Accordingly, the first thought that came to mind was the idea 
of extending the concept of Euler angles to the n-D case. It turned out that, although not well known, Euler 
himself succeeded in using Euler angles to parametrize higher dimensional orthogonal matrices. Euler, 
however, was not concerned with the dynamic case; that is, with the differential equation which describe 
their change in time (neither did he do it for the 3-D case). Lagrange improved Euler’s approach and 
presented it in the first edition of his book on analytic mechanics. We did not adopt this approach because 
of the multitude of the trigonometric functions that one has to trace and compute and because of the com- 
plexity of the differential equations which describe the time change of the extended Euler angles. 

Another popular 3-D parametrization which was considered was the quaternion of rotation. This 
approach did not seem to lead to any solution and was abandoned. The last parametrization which was exa- 
mined was that of Rodrigues. In the vast literature on 3-D methods the elements of this parametrization are 
known as the Gibbs vector or Cayley-Rodrigues parameters; however, the presentation of these parameters 
by Rodrigues in 1840 [7] preceded the work of Cayley who, as a matter of fact, credits Rodrigues with 
their discovery [18, 19]. Rodrigues’ work certainly preceded that of Gibbs who first published his research 
on these parameters in 1884 (see Ref. 9, p. 17). Although it seems that Rodrigues was the first one to 
present them, it turns out, as noted by Jacobi [13] and by Roberson [8], that even these parameters were 
first presented by Euler [20] although in a different form. Ironically, while Rodrigues based his development 
on the, by now, very famous theorem of Euler [21], Euler himself was not aware of the possible use of his 
own theorem in the derivation of these parameters. (The theorem states that any final sequence of 3-D rota- 
tions can be presented by just one rotation about a single fixed axis)*. It is, however, Rodrigues who devel- 
oped the parameters in their present known form. For this reason we refer to their extension to n-D as the 
Extended Rodrigues Parameters. It was shown that the parameters can be conveniently extended to n-D. In 
fact there is nothing that limits their validity to 3-D only. Indeed, the theorems used in the presentation of 
the Extended Rodrigues Parameters in this work do not assume any restriction on the dimensionality of the 
space in which they are used. 

Projecting the 3-D concepts into n-D raises the question of the correct mathematical representation of 
angular rates in spaces whose dimension is not 3. It is shown that angular rate has to be represented by a 
tensor of the second rank, also known as dyadic. The ability to represent angular rate as a vector is unique 
to 3-D. This fact, while known before, was not paid sufficient attention because the vectorial representation 
satisfied the intuition and the practical needs of its users. In other dimensions the vectorial representation 
fails and the use of the dyadic representation is required. Finally it should be pointed out that when, as in 


* As noted by Jacobi, Lagrange too presented this theorem in the first edition of his book on analytic 
mechanics [12] but dropped it from the second edition of this book. 
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our case, matrices are used, the skew-symmetric dyadic which represents angular rate in n-D is simply 
represented by a skew symmetric matrix. 
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